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ABSTRACT 

This  study  focused  on  flow  in  the  unsaturated  or  vadose  zone,  which  forms  a  major  hydrologic 
link  between  the  ground  surface  and  underlying  groundwater  aquifers.  To  properly  understand  its 
role  in  protecting  groundwater  from  surface  and  near  surface  sources  of  contamination,  one  must 
be  able  to  analyze  quantitatively  fluid  flow  in  unsaturated  soils.  The  difficulty  is  that  such  soils 
are  ubiquitously  heterogeneous,  with  hydraulic  properties  that  fluctuate  from  point  to  point  in  a 
seemingly  erratic  manner.  The  common  approach  has  been  to  delineate  this  variation,  and 
analyze  unsaturated  flow  in  randomly  heterogeneous  soils,  deterministically.  Yet  with  increasing 
frequency,  the  popular  deterministic  approach  is  proving  to  be  inadequate.  Our  project  aimed  at 
developing  theoretical  and  computational  methods  to  predict,  in  an  optimum  fashion,  unsaturated 
flow  in  randomly  heterogeneous  soils  under  the  action  of  uncertain  forcing  tenns  and  to  assess 
the  corresponding  prediction  errors.  Previously,  such  predictions  and  assessments  required  the 
conduct  of  numerous  Monte  Carlo  simulations  on  a  fine  grid,  which  was  computationally 
demanding  and  therefore  seldom  used  in  practice.  An  alternative  is  to  conduct  Monte  Carlo 
simulations  on  a  coarse  grid,  which  is  still  computationally  intensive  (due  to  the  need  for  many 
repetitions)  and  leads  to  a  loss  of  accuracy  due  to  the  need  to  average  (upscale)  the  flow 
equations  over  relatively  large  grid  cells.  Our  objective  was  to  avoid  the  need  for  either  Monte 
Carlo  simulation  or  upscaling  by  developing  ways  to  render  predictions  and  uncertainty 
assessments  directly,  in  a  computationally  efficient  and  accurate  manner.  This  final  technical 
report  describes  our  accomplishment  in  the  development  of  two  novel  approaches,  one  based  on 
the  Kirchhoff  transformation  and  the  other  on  a  Gaussian  method  of  approximation.  The  report 
also  describes  some  initial  ideas  about  how  to  extend  the  applicability  of  these  two  approaches  to 
multidimensional  and  transient  flows  in  a  broader  class  of  soils  than  we  have  considered 
previously.  The  report  cites  publications  based  on  research  supported  in  part  by  this  ARO  grant. 

INTRODUCTION 

Our  project  dealt  with  the  effect  of  measuring  randomly  varying  soil  hydraulic  properties  on 
one's  ability  to  predict  unsaturated  flow  subject  to  random  sources  and/or  initial  and  boundary 
conditions.  Our  aim  was  to  develop  theoretical  and  computational  methods  for  the  optimum 
prediction  of  unsaturated  flow  (in  terms  of  pressure  head,  water  content,  flux  and  velocity)  in 
randomly  heterogeneous  soils  under  the  action  of  uncertain  forcing  terms  (boundary  conditions, 
initial  conditions,  sources  and  sinks),  and  assessment  of  the  associated  prediction  errors.  In  the 
past,  such  predictions  and  assessments  required  the  conduct  of  numerous  (hundreds  or 
thousands)  Monte  Carlo  simulations  on  a  fine  computational  grid,  which  is  computationally 
demanding  and  therefore  seldom  used  in  practice.  An  alternative  is  to  conduct  Monte  Carlo 


simulations  on  a  coarse  grid,  which  is  still  computationally  intensive  (due  to  the  need  for  many 
repetitions)  and  leads  to  a  loss  of  accuracy  due  to  the  need  to  average  (upscale)  the  flow 
equations  over  relatively  large  grid  cells.  Our  objective  was  to  avoid  the  need  for  either  Monte 
Carlo  simulation  or  upscaling  by  developing  ways  to  render  predictions  and  uncertainty 
assessments  directly  (in  a  finite  number  of  computational  steps)  in  a  computationally  efficient 
and  accurate  manner.  Whereas  the  Monte  Carlo  method  requires  specifying  the  probability 
distribution  of  soil  parameters  (which,  for  convenience,  are  typically  assumed  to  be  multivariate 
Gaussian  or  log-Gaussian,  mostly  without  direct  evidence),  one  of  the  two  approaches  we  are 
pursuing  (that  based  on  the  Kirchhoff  transformation)  is  free  of  such  distributional  requirements. 

Our  study  focused  on  flow  in  the  unsaturated  or  vadose  zone.  This  zone  forms  a  major 
hydrologic  link  between  the  ground  surface  and  underlying  groundwater  aquifers.  To  properly 
understand  its  role  in  protecting  groundwater  from  surface  and  near  surface  sources  of 
contamination,  one  must  be  able  to  analyze  quantitatively  fluid  flow  in  unsaturated  soils.  The 
difficulty  is  that  such  soils  are  ubiquitously  heterogeneous,  with  hydraulic  properties  (saturated 
conductivity,  porosity,  parameters  of  constitutive  functional  relationships  between  relative 
conductivity,  pressure  head  and  saturation)  that  fluctuate  from  point  to  point  in  a  seemingly 
erratic  manner.  This  erratic  spatial  variability,  together  with  random  errors  of  measurement  and 
interpretation,  renders  the  soil  parameters  uncertain  and  the  corresponding  flow  equations 
stochastic.  In  practice,  the  random  spatial  variability  of  soil  hydraulic  properties,  and  the 
stochastic  nature  of  unsaturated  flow  variables  (pressure  head,  saturation,  flux,  velocity),  are 
often  ignored.  Instead,  the  common  approach  is  to  delineate  the  spatial  variation  of  soil 
properties  deterministically  and  express  the  unsaturated  flow  equations  in  a  similar  fashion. 
Another  popular  assumption  is  that  flow  through  the  vadose  zone  is  largely  vertical  and  one  is 
justified  ignoring  lateral  variations  in  soil  hydraulic  properties  and  flow  variables.  Yet  with 
increasing  frequency,  such  attitudes  are  proving  to  be  counter  productive.  Consider  for 
illustration,  the  case  of  leaking  underground  tanks  at  Hanford. 

A  recent  report  to  Congress  by  the  US  General  Accounting  Office  (GAO/RCED-98-80,  March 
1998)  suggests  that  the  common  practice  of  ignoring  lateral  flow  might  have  misled  the  DOE  to 
believe,  for  many  years,  that  the  vadose  zone  at  Hanford  constitutes  an  effective  barrier  for 
contaminant  migration  between  tank  wastes  and  underlying  groundwater.  The  DOE  had  assumed 
that  wastes  would  move  slowly,  if  at  all,  through  the  vadose  zone,  thereby  obviating  the  need  for 
detailed  studies  of  flow  conditions  in  the  thick  unsaturated  zone  at  Hanford.  The  GAO  report 
cites  evidence  that  may  indicate  otherwise.  Indeed,  in  December  1997  the  DOE  had  announced 
publicly  that  highly  radioactive  wastes  from  previously  leaking  underground  storage  tanks  had 
migrated  all  the  way  down  to  groundwater.  Current  understanding  of  how  wastes  move  through 
the  vadose  zone  to  the  groundwater  has  proven  to  be  inadequate  for  key  technical  decisions  on 
how  to  clean  up  the  wastes  at  Hanford  in  an  environmentally  sound  and  cost-effective  manner. 
To  better  understand  fluid  flow  and  contaminant  transport  processes  in  the  vadose  zone,  one 
must  recognize  that  unsaturated  soils  and  rocks  form  part  of  a  complex  three-dimensional, 
multiphase,  multiscale  heterogeneous  and  anisotropic  hydrogeologic  system.  This  system  does 
not  constitute  a  perfect  sequence  of  horizontal  layers;  if  it  did,  flow  and  transport  rates  would  be 
controlled  by  the  least  permeable  layer  and  would  therefore  be  correspondingly  low.  In  reality, 
unsaturated  medium  properties  vary  spatially  in  a  complex  manner,  which  often  allows  fluids 
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and  contaminants  to  move  around  low-permeability  obstacles  much  faster  than  would  be  possible 
in  the  perfectly  stratified  case.  Preferential  flow  through  high-penneability  channels  and/or  the 
formation  of  unstable  fingers  could  further  enhance  the  rate  of  contaminant  migration  from  a 
source  in  the  vadose  zone  to  the  water  table.  A  panel  of  four  vadose  zone  experts  concluded 
(DOE/RL-97-49,  April  1997)  that  characterization  of  the  vadose  zone  at  Hanford  is  an  essential 
step  toward  understanding  contamination  of  the  groundwater,  assessing  the  resulting  health  risks, 
and  defining  the  concomitant  groundwater  monitoring  program  needed  to  verify  risk 
assessments.  As  flow  and  transport  in  the  heterogeneous  vadose  zone  at  Hanford  are  poorly 
understood,  previous  and  ongoing  computer  modeling  efforts  are  inadequate  and  based  on 
unrealistic  and  sometimes  optimistic  assumptions,  which  render  their  output  unreliable. 

A  recent  study  of  groundwater  and  soil  cleanup  by  the  U.S.  National  Academy  (National 
Research  Council,  1999)  recognizes  that  the  geologic  and  geochemical  characteristics  of  a  site 
have  a  major  influence  on  the  perfonnance  of  subsurface  cleanup  systems.  According  to  this 
study,  the  subsurface  is  usually  highly  heterogeneous  and  characterizing  this  variability  is 
extremely  difficult.  This  heterogeneity  and  difficulty  in  characterization  complicate  the  design  of 
subsurface  cleanup  systems  because  predicting  system  performance  under  such  uncertain 
conditions  is  difficult.  Further,  many  types  of  cleanup  systems,  including  not  only  pump-and- 
treat  systems  but  also  systems  using  in  situ  chemical  oxidation,  biodegradation,  and  other 
processes,  require  the  circulation  of  water,  aqueous  solutions,  or  other  fluids  underground.  The 
physical  heterogeneity  of  the  subsurface  interferes  with  unifonn  delivery  of  fluids  to 
contaminated  locations.  As  a  result,  some  contaminated  zones  will  receive  little  or  no  treatment 
if  a  fluid  is  pumped  in  or  out  of  the  zone.  Technologies  for  treating  DNAPL  source  zones,  and 
dissolved  plumes  emanating  from  DNAPL  sources,  are  limited  primarily  by  geological 
heterogeneities,  which  can  interfere  with  circulation  of  treatment  fluid  and  water  or  can  limit 
access  to  the  subsurface.  Soil  heterogeneity  affects  soil  vapor  extraction  performance  as  air  flows 
most  easily  through  coarse-grained  soils  and  very  little  if  at  all  through  predominantly  clayey 
soils.  Frequently,  volatile  organic  compounds  will  accumulate  preferentially  on  the  surface  of 
and  within  clay  lenses  and  layers,  and  airflow  will  be  minimal  in  the  most  highly  contaminated 
soils.  An  accurate  knowledge  of  geological  heterogeneities  is  vital  for  evaluating  the 
hydrogeological  limits  on  subsurface  contaminant  remediation. 

In  a  January  2000  Editorial  titled  "It's  the  Heterogeneity!”  the  Editor  of  the  most  widely  read 
groundwater  journal  (Wood,  2000)  reminds  his  readers  that  the  heterogeneity  of  chemical, 
biological,  and  flow  conditions  should  be  a  major  concern  in  any  remediation  scenario.  In  his 
view,  many  in  the  groundwater  community  either  failed  to  "get"  the  message  or  were  forced  by 
political  considerations  to  provide  rapid,  untested,  site-specific  active  remediation  technology.  It 
was  their  lack  of  appreciation  of  heterogeneity  that  led  them  to  the  belief  that  they  could 
remediate  aquifers  by  simply  pumping  out  and  treating  offending  solute.  "It's  the  heterogeneity," 
and  it  is  the  Editor's  guess  that  the  natural  system  is  so  complex  that  it  will  be  many  years  before 
one  can  effectively  deal  with  heterogeneity  on  societally  important  scales. 

The  purpose  of  our  work  under  this  ARO  grant  was  to  help  accelerate  the  process.  Though  the 
complex  and  uncertain  nature  of  subsurface  flow  conditions  is  now  widely  recognized,  there 
does  not  yet  appear  to  be  a  satisfactory  way  to  characterize  and  quantify  them  mathematically 
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and  computationally.  Our  goal  was  to  develop  a  mathematical  framework,  and  computational 
algorithms  that  help  materially  advance  the  corresponding  state  of  science  and  technology.  We 
pursued  this  goal  through  (a)  the  development  of  stochastically-derived  deterministic 
"conditional  moment  equations"  that  allow  optimum  unbiased  prediction  of  flow  in  randomly 
heterogeneous  unsaturated  soils  on  a  multiplicity  of  spatial  scales,  under  the  action  of  uncertain 
forcing  terms,  as  well  as  assessment  of  the  corresponding  prediction  errors,  and  (b)  the 
development  of  associated  analytical  and  computational  methods  of  solution.  Rather  than 
requiring  repeated  Monte  Carlo  simulations  on  a  fine  grid,  our  computational  methods  are 
designed  to  yield  a  solution  in  one  single  simulation  on  a  relatively  coarse  grid.  Rather  than 
working  with  upscaled  quantities  which  are  often  difficult  to  justify  theoretically,  or  compare 
with  measurements,  all  quantities  that  enter  into  our  equations  are  defined  on  a  scale  compatible 
with  potentially  available  field  data  (their  support  scale). 

ACCOMPLISHMENTS 

The  grant  funded  a  doctoral  student,  Mr.  Donghai  Wang,  who  plans  to  complete  his  Ph.D. 
dissertation  in  spring  2003.  The  grant  also  provided  partial  funds  for  the  Co-PI,  Dr.  Daniel  M. 
Tartakovsky,  of  the  Scientific  Computing  Group  within  the  Theoretical  Division  at  Los  Alamos 
National  Laboratory.  In  addition  to  the  Principal  Investigator,  persons  not  funded  by  the  grant  in 
any  major  way,  who  however  have  collaborated  with  us,  include  Dr.  Oma  Amir  and  Dr.  Zhiming 
Lu.  The  latter  two  have  been  funded  primarily  by  a  National  Science  Foundation  grant  that 
expired  on  January  31,  2001.  Dr.  Amir  has  completed  her  doctorate  in  December  1999  and  Dr. 
Lu  in  May  2000.  Both  the  PI  and  Dr.  Tartakovsky  have  served  as  members  of  Dr.  Amir’s  and  Dr. 
Lu’s  doctoral  committees.  We  describe  briefly  the  research  and  relevant  work  products 
completed  by  members  of  this  team  to  date. 

Laying  the  Groundwork 

Our  research  was  founded  on  two  papers  that  have  appeared  in  the  fall  of  1999.  In  the  paper  by 
Neuman  et  al.  (1999)  we  lay  the  foundations  of  a  new  conditional  moment  approach  for  the 
solution  of  stochastic  unsaturated  flow  equations  with  random  parameters  and  uncertain  forcing 
tenns.  The  paper  shows  that  the  approach  leads  to  detenninistic  equations  for  the  conditional 
moments,  which  can  be  solved  by  standard  numerical  methods,  thereby  obviating  the  need  for 
Monte  Carlo  simulations.  The  parameters  in  these  equations  may  be  local  (depending  on  one 
point  in  space-time)  or  nonlocal  (depending  on  two  such  points).  As  they  are  conditional  on 
measurements,  the  parameters  are  not  unique  properties  of  the  soil  but  vary  with  the  underlying 
database.  The  conditional  mean  solution  constitutes  an  optimum  unbiased  predictor  of  the 
otherwise  unknown  state  of  the  system,  and  conditional  second  moments  provide  a  measure  of 
the  corresponding  prediction  uncertainty.  Since  all  moments  are  defined  on  the  same  consistent 
measurement  (support)  scale  co  as  the  data,  there  is  no  need  for  upscaling,  though  one  can  easily 
integrate  the  conditional  mean  solution  in  space-time,  if  one  so  desires.  As  the  conditional  mean 
solution  is  smooth  relative  to  its  random  counterpart,  it  can  in  principle  be  resolved  on  a 
numerical  grid  which  is  coarser  than  that  typically  required  for  the  Monte  Carlo  simulation  of 
random  fields. 
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The  paper  of  Neuman  et  al.  (1999)  describes  two  methods  for  the  development  of  conditional 
moment  equations  and  their  solution.  One  method  is  based  on  the  Kirchhoff  transfonnation  and 
the  other  on  a  Gaussian  closure  approximation.  The  method  using  Kirchhoff  transformation  is 
described  more  fully  in  a  paper  by  Tartakovskv  et  al.  (1999).  There  we  demonstrate  that  the 
Kirchhoff  transformation  fully  linearizes  the  stochastic  steady  state  unsaturated  flow  equation  in 
the  absence  of  gravity,  and  does  the  same  for  flow  with  gravity  when  the  Gardner  model,  K  =  Ks 
exp  {ay/),  applies;  here  K  is  unsaturated  hydraulic  conductivity,  Ks  is  saturated  hydraulic 
conductivity,  or  is  a  positive  exponent,  and  y/  is  (negative)  pressure  head.  Consequently,  the 
method  yields  exact  conditional  mean  flow  equations  for  both  situations,  as  well  as  exact 
equations  for  the  conditional  variance-covariance  of  pressure  head  and  flux.  Both  sets  of 
equations  are  nonlocal  (integro-differential)  in  that  they  contain  parameters  depending  on  more 
than  one  point  in  space.  Both  the  local  (depending  on  one  point  in  space)  and  nonlocal 
parameters  in  our  moment  equations  are  conditional  on  data  and  therefore  nonunique.  The 
conditional  mean  solution  constitutes  an  optimum  unbiased  predictor  of  the  otherwise  unknown 
state  of  the  system,  and  conditional  second  moments  provide  a  measure  of  the  corresponding 
prediction  uncertainty. 

The  linear  Kirchhoff-transfonned  stochastic  flow  equations  yield  exact  conditional  moment 
equations  which,  however,  cannot  be  solved  without  a  closure  approximation.  The  closure 
approximation  we  use  is  based  on  perturbation  analysis.  As  such,  it  is  nominally  limited  either  to 
mildly  heterogeneous  soils,  or  to  strongly  heterogeneous  soils  in  which  hydraulic  properties  have 
been  measured  with  sufficient  accuracy,  at  a  sufficiently  large  number  of  points,  to  allow 
estimating  them  everywhere  else  in  the  soil  with  a  relatively  low  degree  of  uncertainty.  We  shall 
see  later  that,  in  reality,  our  perturbation  approach  works  well  even  in  strongly  heterogeneous 
soils  in  the  absence  of  such  conditioning. 

The  paper  of  Tartakovskv  et  al.  (1999)  demonstrates  rigorously  that  the  concept  of  effective 
hydraulic  conductivity  does  not  generally  apply  to  statistically  averaged  unsaturated  flow 
equations  except  when  they  are  unconditional  and  flow  is  driven  solely  by  gravity.  It  points  out 
that  all  conditional  parameters  and  moments  in  our  equations  are  smooth  relative  to  their  random 
counterpart  and  can  therefore  be  resolved,  in  principle,  on  a  numerical  grid  which  is  coarser  than 
that  typically  required  for  the  Monte  Carlo  simulation  of  random  fields.  The  paper  proceeds  to 
develop  analytical  solutions  for  the  Kirchhoff  potential,  pressure  head  and  their  variances  under 
vertical  infiltration,  without  conditioning,  to  second  order  of  approximation  in  the  standard 
deviation  (first  order  in  the  variance)  of  natural  log  saturated  hydraulic  conductivity.  It  then 
compares  these  with  Monte  Carlo  results  obtained  by  solving  the  stochastic  Richards  equation 
numerically.  Our  second  order  approximations  are  generally  far  superior  to  zero  order 
approximations,  and  the  variance  of  pressure  heads  compares  much  better  with  Monte  Carlo 
values  than  does  the  variance  of  Kirchhoff  potentials.  Both  the  analytical  pressure  head  and  its 
variance  compare  well  with  Monte  Carlo  results  for  natural  log  conductivity  variances  at  least  as 
large  as  1 .  This  accords  well  with  theoretical  analysis  (presented  in  the  paper)  which  shows  that 
our  analytical  solution  remains  asymptotic  for  input  variances  as  large  as  2. 
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The  Gaussian  method  of  approximation,  introduced  by  Neuman  et  al.  (1999),  is  more  general 
than  the  method  based  on  Kirchhoff  transfonnation  in  that  it  imposes  fewer  restrictions  on  the 
functional  form  of  constitutive  relationships  between  unsaturated  hydraulic  conductivity, 
pressure  head  and  saturation.  It  however  requires  assuming  that  the  reference  pressure  head  is 
multivariate  Gaussian  (or  log-Gaussian)  about  its  conditional  mean. 

The  Kirchhoff  method  of  solution  has  been  developed  and  explored  further  by  Drs.  Lu  and 
Tartakovsky.  The  Gaussian  method  has  been  the  research  focus  of  Dr.  Amir  and  Mr.  Wang.  We 
present  a  brief  description  of  their  accomplishments  and  findings. 

Numerical  Analysis  Based  on  Kirchhoff  Transformation 

In  May  2000,  Zhiming  Lu  completed  his  doctoral  dissertation  (Lu,  2000)  on  nonlocal  finite 
element  analysis  of  conditional  steady  state  unsaturated  flow  in  bounded,  randomly 
heterogeneous  soils  using  the  Kirchhoff  transformation.  The  highlights  of  his  work  are  described 
in  a  paper  by  Lu  et  al.  (2002)  published  earlier  this  year  in  the  archival  journal  Water  Resources 
Research. 

Dr.  Lu’s  work  consists  of  the  development,  and  computer  implementation,  of  a  finite  element 
algorithm  for  the  prediction  of  steady  state  unsaturated  flow  in  the  vertical  plane.  He  considers 
flow  in  a  bounded,  randomly  heterogeneous  soil  profile  under  the  influence  of  random  forcing 
terms.  His  aim  is  to  predict  pressure  head  and  flux  at  each  point  in  the  two-dimensional  vertical 
profile  without  resorting  to  Monte  Carlo  simulation,  upscaling  or  linearization  of  the  constitutive 
relationship  between  unsaturated  hydraulic  conductivity  and  pressure  head.  To  achieve  this,  he 
represents  the  latter  relationship  through  Gardner's  exponential  model  (described  earlier), 
treating  its  exponent  a  as  a  random  constant  and  saturated  hydraulic  conductivity,  Ks,  as  a 
spatially  correlated  random  field.  This  allows  him  to  linearize  the  steady  state  unsaturated  flow 
equations  by  means  of  the  Kirchhoff  transformation,  integrate  them  in  probability  space,  and 
obtain  exact  integro-differential  equations  for  the  conditional  mean  and  variance-covariance  of 
transformed  pressure  head  and  flux,  in  the  manner  of  Tartakovsky  et  al.  (1999).  Expansion  of  the 
nonlocal  conditional  moment  equations  in  powers  of  crY  and  oa ,  which  represent  measures  of 
the  standard  estimation  errors  of  saturated  natural  log  hydraulic  conductivity  Y=  In  Ks  and  J3=  In 
a,  respectively,  leads  to  a  set  of  recursive  closure  approximations.  Zhiming  solves  these 
approximate  moment  equations  by  finite  elements,  to  second-order  of  approximation,  for 
superimposed  mean  uniform  and  divergent  flow  regimes. 

As  the  conditional  mean  quantities  are  generally  smoother  than  their  random  counterparts,  the 
recursive  moment  equations  can  be  solved  (in  principle)  on  a  relatively  coarse  grid  without 
upscaling.  Dr.  Lu,  however,  uses  a  fine  grid  to  compare  his  nonlocal  finite  element  solution  with 
conditional  and  unconditional  Monte  Carlo  simulations,  conducted  on  the  same  grid  by  standard 
finite  elements.  Dr.  Lu’s  comparison  demonstrates  that  his  direct  finite  element  solution  of  the 
moment  equations  is  highly  accurate  for  mildly  heterogeneous  soils  and  works  well  for  soils  that 
are  moderately  to  strongly  heterogeneous. 
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Extension  of  Approach  Based  on  Kirchhoff  Transformation  to  Soils  with  Spatially  Varying 
Exponent  a 

Our  analytical  and  numerical  approaches  based  on  the  Kirchhoff  transfonnation  treat  the 
exponent  a  in  the  Gardner  (1958)  constitutive  relation  K  =  Ks  exp  {ay/)  as  a  random  constant 
that  does  not  vary  in  space.  In  a  recent  Journal  of  Hydrology  article  ( Tartakovskv  et  a/.,  2002)  we 
describe  a  way  to  relax  this  requirement  by  allowing  a  to  be  a  statistically  homogeneous  random 
field.  Our  approach  utilizes  the  “partial  mean-field”  concept,  according  to  which  the  random 
field  O'  ( v)  is  replaced  by  its  spatially  varying  conditional  ensemble  mean  while  saturated 

hydraulic  conductivity  remains  a  random  field.  The  Kirchhoff  transformation  can  then  be  applied 
to  the  resulting  (nonlinear)  stochastic  partial  differential  equation  in  a  manner  similar  to  that  of 
our  earlier  analyses.  The  accuracy  of  this  approach  depends  on  a  complex  interplay  between  the 
statistical  parameters  of  O')*)  (mean,  variance  and  correlation  scale),  an  issue  explored  in  the 
above  paper. 

Analysis  Based  on  Gaussian  Approximation 

In  December  1999  Orna  Amir  completed  her  doctoral  dissertation,  titled  "Gaussian  Analysis  of 
Unsaturated  Flow  in  Randomly  Heterogeneous  Porous  Media"  (Amir,  1999).  Based  on  the 
assumption  that  pressure  head  y/  is  multivariate  Gaussian  about  its  conditional  mean,  Dr.  Amir 
was  able  to  derive  governing  equations  for  the  mean  and  variance  of  y/  without  linearizing  either 
the  constitutive  relation  between  unsaturated  hydraulic  conductivity  and  pressure  head,  or  the 
governing  unsaturated  flow  equations  (so  that  the  governing  moment  equations  remain  nonlinear, 
as  is  the  underlying  stochastic  Richards’  equation).  Contrary  to  all  other  known  solutions  of  the 
stochastic  unsaturated  flow  problem,  our  Gaussian  approach  places  no  obvious  restrictions  on  the 
variance  of  the  corresponding  constitutive  parameters.  This  is  evident  from  Dr.  Amir’s 
computational  results. 

Dr.  Amir  illustrated  the  application  and  effectiveness  of  the  Gaussian  closure  approximation  by 
developing  a  closed  system  of  coupled  nonlinear,  ordinary  differential  equations  for  the  first  and 
second  moments  of  pressure  head  under  one-dimensional  steady  state  unsaturated  flow  through  a 
randomly  stratified  soil.  Her  equations  are  written  (by  choice,  not  necessity)  for  unsaturated 
hydraulic  conductivity  that  varies  exponentially  with  pressure  head,  where  now  the  exponent  a 
is  not  a  random  constant  (as  it  was  in  the  case  of  the  Kirchhoff  transformation)  but  a  spatially 
varying  random  field.  Rather  than  treating  the  soil  as  a  continuum,  Dr.  Amir  found  it  helpful  to 
represent  it  by  a  discrete  assembly  of  layers,  each  of  which  has  uniform  but  random  properties  Y 
=  In  Ks  and  fd  =  In  a ,  which  however  are  auto-  and  cross-correlated  between  the  layers. 

Dr.  Amir  solved  her  one-dimensional,  steady  state  nonlinear  moment  equations  numerically  and 
compared  the  results  with  those  obtained  by  Monte  Carlo  simulation,  based  on  an  existing 
analytical  solution  of  Richards'  equation  for  this  case.  The  comparison  shows  excellent 
agreement  between  the  two  sets  of  results  over  a  remarkably  broad  range  of  constitutive 
parameters.  A  paper  that  describes  this  part  of  Dr.  Amir’s  work  has  appeared  in  the  journal 
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Transport  in  Porous  Media  ( Amir  and  Neuman ,  2001a).  An  earlier,  less  complete  version  has 
appeared  in  a  book  published  as  a  Geological  Society  of  America  Special  Paper  ( Amir  and 
Neuman,  2000). 

Extension  of  Gaussian  Approach  to  Transient  Flow 

In  her  dissertation,  Dr.  Amir  extended  her  one-dimensional  Gaussian  solution  to  the  case  of 
transient  flow  with  an  exponential  constitutive  relationship  between  saturation  and  pressure  head. 
Upon  comparing  her  mean  solution  with  that  obtained  by  the  Monte  Carlo  method  (this  time 
through  numerical  solution  of  the  one-dimensional  Richards'  equation),  Dr.  Amir  found  that  the 
two  agree  very  well  for  a  wide  range  of  constitutive  parameters.  However,  she  found  a  less 
satisfactory  agreement  between  the  variances  and  covariances  of  pressure  head  obtained  by  the 
two  methods. 

Though  Dr.  Amir  is  now  in  Israel,  she  continues  to  collaborate  with  us  on  this  project.  During  the 
last  year,  Dr.  Amir  was  able  to  improve  substantially  the  quality  of  pressure  head  variance  and 
covariance  assessments  based  on  her  one-dimensional  Gaussian  approach.  A  paper  summarizing 
her  most  recent  work  on  this  topic  has  been  submitted  for  publication  in  the  journal  Transport  in 
Porous  Media  (Amir  and  Neuman,  2002b). 

Extension  of  Gaussian  Approach  to  Multidimensional  Flow 

To  lay  the  groundwork  for  a  multi-dimensional  application  of  our  Gaussian  approach,  Amir 
(1999)  proposed  solving  steady  state  flow  in  a  two-dimensional  domain  by  finite  elements  by 
representing  the  soil  as  a  checkerboard  of  square  elements,  each  having  uniform  random 
constitutive  properties  Y  =  In  Ks  and  fd  =  In  a ,  which  however  are  auto-  and  cross-correlated 
between  the  elements.  The  computational  implementation  of  this  finite  element  algorithm  has 
become  the  domain  of  a  doctoral  student  supported  by  this  grant,  Mr.  Donghai  Wang. 

Donghai  has  formulated,  developed  and  implemented  a  finite  element  algorithm  based  on  the 
above  idea.  He  has  applied  his  algorithm  to  two-dimensional  flow  in  a  bounded  vertical  domain 
under  coupled  mean  uniform  and  convergent  flows,  and  compared  his  results  with  those  of 
standard  Monte  Carlo  simulations.  His  excellent  results  are  summarized  in  a  paper  presented  at 
(and  included  in  the  Proceedings  of)  the  Fourteenth  International  Conference  on  Computational 
Methods  in  Water  Resources  that  took  place  in  Delft,  The  Netherlands,  in  June  2002  (Wans  et 
al,  2002).  Mr.  Wang  is  presently  completing  his  doctoral  dissertation,  which  he  plans  to  defend 
in  spring  2003.  We  plan  to  prepare  one  or  more  journal  articles  based  on  his  dissertation  in 
spring  and  summer  2003. 

Following  is  a  brief  technical  description  of  Mr.  Wang’s  work. 

Consider  steady  state  unsaturated  flow  in  a  bounded  domain  Q  governed  by  mass  continuity  and 
Darcy’s  law, 
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-V-q(x)  +  /(x)  =  0  q(x)  =  -A(x,^)V(^  +  x3) 

subject  to  boundary  conditions 

if/(x)  =  vF(x)  on  Tfl 

-q(x)-n(x)  =  0(x)  on 


(1) 


(2) 


Here  V  is  gradient  operator  with  respect  to  the  spatial  position  vector  x ,  q  is  flux,  /  is  a  source 
tenn,  K  is  a  spatially  correlated  random  hydraulic  conductivity  field,  i/7  is  pressure  head,  x3  is 
the  vertical,  To  is  Dirichlet  boundary,  r  v  is  Neumann  boundary,  and  n  is  a  unit  outer  normal  to 
the  boundary.  The  forcing  terms/  'F  ,  Q  are  random  and  mutually  uncorrelated.  All  quantities 
are  defined,  and  measurable,  on  a  bulk  support  volume  co  that  is  small  compared  to  the  flow 
domain  Q  .  Flow  takes  place  under  strictly  unsaturated  conditions  such  that  i//  <  0.  The  random 
nature  of  K  and  the  forcing  tenns  render  (1)  -  (2)  stochastic. 

We  represent  hydraulic  conductivity  using  the  Gardner  (1958)  constitutive  model 

K(x,i//)  =  Ks(x)Kr(x,i//)  Kr(x,i//)  =  ea(x)¥(x)  (3) 


where  Ks  is  saturated  hydraulic  conductivity,  Kr  is  relative  conductivity,  and  a  is  a  positive 
exponent.  The  flux  in  (1)  and  (2)  then  becomes 

q  {x)  =  -Ksea*V{y  +  x,)  (4) 


We  treat  7(x)  =  In  Ks  (x)  as  a  correlated  Gaussian  random  field  and  A  =  -  In  a  as  a  normally 
distributed  random  variable.  Using  (4),  we  rewrite  (1)  -  (2)  in  dimensionless  form  as 


0  = 


e¥(x)S7\jxy/(x)  +  (a) 


(5) 


^'(x)  =  vF(x)  onT^ 

_^(^)ed-v>v[Q,^(x)  +  (Q,)x3].n(x)  =  g(x)  on/ 

{a) 

where  (a)  is  ensemble  mean  of  a ,  V  is  gradient  operator  with  respect  to  dimensionless  x,  and 
y/  and  x  are  dimensionless  variables  defined  as 

yz-y/la  x  =  x/(a)  (7) 


9 


Considering  that(er)  =  e  '/('+CT,/2  where  (d)  is  the  mean  of  A  and  a]  is  its  variance,  we  can  rewrite 
(5)  -  (6)  as 


_ 


0  =  e 


V- 


Y+{A)-Aj  2 


evV\e  y/  +  e 


^(jc)  =  'P(jc)  on  rD 

_ei+C)-^/2evy  ^e~Ay/  +  e~^<TAl2xi  j-n(x)  =  Q(x)  on 


(8) 

(9) 


We  take  the  flow  domain  to  be  a  checkerboard  of  R  densely  spaced,  nonoverlapping  subdomains 
Q  r,r=  1,2,  . .  .R,  within  each  of  which  Y  and  A  are  random  constants,  Yr  and  Ar: 

f(x)  =  {^,^etl,.Vr}  ^(x)  =  {4.,*eQ,.Vr}  (10) 


Multiplying  (8)  by  a  detenninistic  weight  function  fa ,  integrating  over  the  global  domain  Q , 
rewriting  as  the  sum  of  integrals  over  R  sub-domains  Q;. ,  and  taking  ensemble  mean  yields 


0  -  £ e{Ar)~aX  n  f  V  •  leY+{A)~aX  newV  (e~Ayr  +  e~{A)+<  12 x3  )\  fa,  AD. 

r= 1  n,  '  '  A 

Applying  Green’s  identity,  then  rewriting  yr  as  (iff)  +  (//  and  Y as  (Y)  +  Y' ,  (1 1)  becomes 


(ID 


r= 1 

+  ^e(Yr)+(Ar)-Ar 


J  le'l'  +Y'y[e  A'  (yr)  +  e  A  y/'  +  e*'  nxi  j 
J  fae^  L ev+Y  v(e~A  (i/r)  +  e~A  y/'  +  eaAr  '2 

rnnr  '  V 


•Vfa  dYi, 


(12) 


Treating  yr,  Y  and  A  as  Gaussian  and  defining 
p  =  e(y)  (f'+Y'-A')  £  =  (yr'A,)  +  (YA' ) 

allows  simplifying  (12)  as 

0  =  - a,r  | (Vp  +  pe^i^-Vfa  dCl  +  y'e^'1^'1  aAr  J  fa(s/p+  peH^ndT 

r= i  n,  r= i  rnn, 

where  i3  =  (0,0,  l)7  .  Approximating  p  and  ^by 


(13) 


(14) 
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(15) 


P'^pj..  C  =  t.CJ. 


where  pm  and  C,m  are  values  of  p  and  C,  at  nodes  m  =  1,  2,  ...  N  and  now  serve  as  Lagrange 
interpolation  functions,  then  using  (9),  yields  the  Galerkin  finite  element  scheme 


°=-z 


R  N 

(Y,H  Ar)-o\ 


r-\ 

R 


«=1  nr 


Z  Pm  J  I  V  K  ■  S7</>„  +  1  e{K)+M~aA'  Z  Pm  J  </>y  </>,„  ■  n  dT 


»=i  rDnn,. 


(Yr)+(Ar 

e 

r=t  rflnn 


+z* 


R 


'■  f  J  <j>nQdY 


M- 

\ 

r=l 


(16) 

where  P  is  p  at  the  r n  fltf. .  Setting  n  =  1,2,  . . .  N  yields  a  system  of  N  equations  in  N  unknown 
values  of pn.  However,  the  equations  are  nonlinear  due  to  their  dependence  on  £ 


To  derive  a  complementary  system  of  equations  for  we  multiply  (1 1)  by  A’  and  (j)n,  integrate 
over  Q  ,  rewrite  as  a  sum  of  integrals  over  Llr ,  take  ensemble  mean,  apply  Green’s  identity  and 
account  for  (9)  to  obtain 


o  =  -Z 


r= 1 

R 


+z 


jYrHA] 


(Yr)+(Ar)~< 


Z  p£i  1 v (M- yvt'dci- a2 Z  p,  j  V 6  •  v h  dn  +  X  J 


/3  •  V  (f)n  dQ. 


i,j= 1  Q,. 

N 


i= 1  Q,. 

N 


*.7=1  fi,. 


Z  j  </>„’ v  {Mj )  •  n  rfr  -o-j  Z  Pi  j  ‘t’n  V  </>,  •  n  c/r 

\yj=i  rflnn,  1=1  rflnn,  y 


+  Ze<y'>+^)-<  |  ^((^V^e^Vnrfr-Xe41'*  j  </>,,  {A’Q)dY 

r=1  rflnn,  r=1  rKnn, 


(17) 

Setting  n  —  1,2,  ...TV  yields  a  system  of  Af  equations  in  Af unknown  values  of  •  The  coupled 
nonlinear  equations  (16)  -  (17)  are  solved  simultaneously  by  (in  our  case  Picard)  iteration  for  p 
and  £ 


To  compute  the  covariance  function  Cv  (jc, x)  =  (y' {x)y/' {xfj  of  dimensionless  pressure  head 
(x  and  x  being  two  arbitrary  points  in  space),  we  first  require  an  equation  for  the  mixed 
moment  .  Multiplying  (8)  by  Y'  and  following  the  same  procedure  as 

before  yields  a  system  of  A^  linear  equations  for  nodal  values  of  <p. 


o  =  -Z 


Pr)+{A)-0-\ 


Z <Pm  ~ J  ( <t>yp  +  pv <t>m  +  p</>,„ech )-^</>„dn+  j  <t>n (</>yP + Pv <f \n 


rnn. 


\ 

+  p<j>meci3 )  ■  n  dT 

7 


+  yg(WW 


( 


-J[(^-(ra'»Vp  +  <T^efi3]-V^  dn+  j  <A[(cry-(7fi'))Vp 

V.  n,  rnnr 


\ 

( Type 41  i3 1  ■  n  dT 


m 


ii 


Multiplying  (8)  by  i//'  =  i//'(x)  and  following  a  similar  procedure  yields  a  system  of  NxN 
equations  for  Cmm  =  C,;/  (xn ,  x- ) , 


°  =  -Ze<1'>+(4><T''ZC™  “J  {</,mVP  +  PV</>m+P</>„/h)-V</>n  dn+  J  <t>n{</>yp  +  pv  </>,,, +P<l>me‘:h)-n  dr 


r= 1 
R 

-5> 

r=l 


rnn. 


-\\y{p(*P-ili)+VPe(i 3  'Y<j)„dQ.+  j  [v (/?(<£>- )  +  <£> peH^  -nt/r 


rnn. 


(19) 

Computation  is  facilitated  by  the  fact  that  (18)  and  (19)  contain  identical  coefficient  matrices  and 
the  symmetric  nature  of  the  matrix  Cm  ~ . 


Additional  details  about  these  derivations  can  be  found  in  Amir  (1999). 

The  variance  cr  is  obtained  from  the  covariance  C://  (x,x)  by  setting  x  =  x  .  Finally,  the  mean 
solution  is  obtained  via 


(ip)  =  In  p-(<jy  /  2  +  crj/2  +  cC/2  +  ^>) 


(20) 


For  illustration  purposes  we  consider  flow  in  a  vertical  plane  of  size  4x8  (all  terms  being  given  in 
arbitrary  consistent  units)  having  impermeable  side  boundaries  (Figure  1).  A  constant  deterministic 
flux  Qb  =  0.01  is  prescribed  at  the  top  and  zero  pressure  head  at  the  bottom.  A  point  source  of 
magnitude  Qs  =  1  causes  the  otherwise  near-unifonn  mean  flow  to  become  locally  divergent  in  the 
domain  interior.  For  simplicity,  Y  and  A  are  taken  to  be  spatially  and  mutually  uncorrelated.  We 
solve  the  problem  using  our  finite  element  Gaussian  closure  algorithm,  and  (for  comparison)  by 
5,000  Monte  Carlo  simulations  using  standard  finite  elements,  on  a  grid  of  20  x  40  square  elements 
with  bilinear  weight  and  interpolation  functions. 
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Parameters: 

a =1.0 
a =0.01 

NMC  =  5000 
Ax,  =  0.2 
Ax2  =0.2 

■  Point  source 


Figure  1.  Problem  definition  and  associated  grid 


Figure  2  depicts  contours  and  profiles  of  mean  dimensionless  pressure  head  and  its  variance 
obtained  by  the  two  methods  for  a  medium  having  uniform  but  random  Y  and  A  values 
characterized  by  (Y\  =3 ,  <j\  =  2,  (A)  =  0,  and  crj  =0.02.  The  two  solutions  are  seen  to  agree 

very  well  despite  the  relatively  large  value  of  <j\ . 

Figure  3  shows  how  the  mean  and  variance  of  Y  and  A  vary  spatially  in  two  cases  for  which  we 
present  solutions  in  Figures  4  and  5.  In  both  cases,  the  Gaussian  Closure  and  Monte  Carlo 
methods  yield  virtually  identical  dimensionless  mean  pressure  head  values  but  slightly  different 
variances.  Overall,  we  consider  these  results  to  be  very  good. 

Our  solution  of  (16)  -  (17)  converges  in  3  to  4  iterations.  This,  and  the  fact  that  we  need  to  solve 
our  Gaussian  closure  equations  only  once,  helps  explain  why  our  solution  has  taken  only  about 
one  fourth  the  time  required  for  the  completion  of  5,000  standard  Monte  Carlo  simulations 
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Figure  2.  Contours  and  profiles  of  mean  pressure  head  (A,  B)  and  variance  of  pressure  head  (C, 
D)  obtained  by  Monte  Carlo  (MC,  solid)  and  Gaussian  closure  (GC,  dash-dot)  for  homogeneous 
domain  with  <Y>  =  3,  a*  =  2  and  <A>  =  0,  crj  =  0.  02 


Figure  3.  Variation  of(f)  and  crj  in  cases  of  Figure  4  (A  and  B)  and  Figure  5  (C  and  D).  In  both 

cases  al  =  1  and  (^4)  =  0. 
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Figure  4.  Contours  and  profiles  of  dimensionless  mean  pressure  head  (A-B)  and  variance  of 
dimensionless  pressure  head  (C-D)  obtained  by  Monte  Carlo  (MC,  solid)  and  Gaussian  closure 
(GC,  dash-dot)  for  parameters  defined  in  Figure  3A-B. 


Figure  5.  Contours  and  profiles  of  dimensionless  mean  pressure  head  (A-B)  and  variance  of 
dimensionless  pressure  head  (C-D)  obtained  by  Monte  Carlo  (MC,  solid)  and  Gaussian  closure 
(GC,  dash-dot)  for  parameters  defined  in  Figure  3C-D. 

Advances  in  Modeling  Soil  Constitutive  Relations 

A  paper  by  Assouline  and  Tartakovskv  (2001)  on  this  topic  has  recently  been  published  in  Water 
Resources  Research.  It  describes  the  development  of  a  new  two-parameter  expression  for 
relative  hydraulic  conductivity  of  partially  saturated  soils.  The  new  expression  is  based  on  a 
premise  by  Assouline  et  al.  (1998)  that  the  probability  distribution  of  particle  volumes  in  natural 
soils  results  from  a  series  of  sequential  fragmentations.  These  fragmentations  are  caused  by 
cyclical  wetting  and  drying;  physical,  chemical  and  biological  processes;  and  cultivation 
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practices.  Assouline  assumes  the  fragmentation  process  to  be  uniform  and  random,  and  the 
probability  of  particle  fragmentation  to  be  proportional  to  its  volume.  This  leads  to  a  probability 
distribution  of  soil  particles  that  is  asymptotically  exponential.  The  latter  is  a  particular  case  of 
the  Weibull  distribution.  Assouline  and  Rouault  (1997)  and  Rouault  and  Assouline  (1998)  have 
established  a  power  relationship  between  particle  volume  and  pore  volume.  It  implies  that  pore 
volume  distribution  is  described  by  the  general  Weibull  model.  Coupling  this  with  the  capillary 
law  has  allowed  Assouline  et  al.  (1998)  to  express  the  relationship  between  effective  saturation 
and  pressure  head  in  terms  of  pressure  head  at  the  wilting  point  and  two  empirical  parameters, 
which  are  determined  by  fitting  the  function  to  experimental  data.  Assouline  and  Tartakovsky 
extended  the  approach  by  deriving  a  corresponding  relationship  between  relative  hydraulic 
conductivity,  effective  saturation  and  pressure  head.  Upon  fitting  their  relative  conductivity 
model  to  data  representing  various  soil  types,  the  authors  found  that  it  fits  these  data  better  than 
do  the  widely  used  models  of  Brooks  and  Corey  (1964)  and  van  Genuchten  (1980). 

SUMMARY  OF  ACCOMPLISHMENTS 

The  goal  of  this  work  was  to  advance,  as  far  as  possible,  our  ability  to  render  reliable  predictions 
of  flow  in  randomly  heterogeneous  soils  under  conditions  of  uncertainty  in  a  computationally 
efficient  manner,  and  to  assess  the  uncertainty  of  these  predictions.  We  have  developed  two 
major  techniques  to  accomplish  this  goal,  one  based  on  the  Kirchhoff  transformation  and  the 
other  on  a  Gaussian  closure  approximation.  Both  techniques  presently  utilize  the  finite  element 
method  to  solve  two-dimensional  steady  state  unsaturated  flow  problems  with  gravity  in  the 
presence  of  arbitrary  source  and  boundary  terms.  We  have  started  working  on  methods  to  extend 
these  techniques  to  transient  flows  in  soils  having  arbitrary  constitutive  properties. 
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